Description

Adds support for cluster-robust (sandwich) variance estimation in both est_seroincidence() and est_seroincidence_by() to account for within-cluster correlation in clustered sampling designs (e.g., household, school-based surveys). Includes support for multi-level clustering (e.g., schools nested within districts).

Changes Made

New Parameters in est_seroincidence() and est_seroincidence_by()

• cluster_var: Cluster identifier variable name(s). Can be a single variable (character string) or multiple variables for multi-level clustering (e.g., c("school", "classroom"))
• stratum_var: Stratum identifier variable name (optional)
• sampling_weights: Reserved for future implementation

Variance Calculation

• Added .compute_cluster_robust_var() implementing sandwich estimator (V = H⁻¹BH⁻¹)
• Computes cluster-level scores via numerical differentiation
• summary.seroincidence() automatically uses cluster-robust variance when clustering detected
• Fixed issue where cluster-robust SEs caused [] notation in column names (SE[,1] instead of SE)
• Added se_type column to summary output indicating "standard" or "cluster-robust"
• Supports multi-level clustering by accepting multiple cluster variables

Code Organization

• Created .validate_cluster_params() helper function to extract cluster and stratum validation logic from main function
• Improved code maintainability by decomposing validation checks into separate subfunction
• Added version management guidelines to .github/copilot-instructions.md documenting requirement to keep dev version one past main branch
• Documented vignette subfile policy in .github/copilot-instructions.md explaining header placement rules and _ prefix naming convention

Bug Fixes

• Fixed est_seroincidence_by() to properly pass cluster and stratum variables through stratified analyses
• Updated stratify_data() to preserve cluster/stratum columns during data stratification (previously these columns were dropped, causing errors when using clustering with est_seroincidence_by())
• Fixed column naming issue where cluster-robust standard errors caused [] notation in column names
• Fixed vignette rendering error in clustering examples by ensuring noise_params are filtered to match pop_data country (Pakistan)

Tests and Documentation

• Comprehensive test suite in test-cluster_robust_se.R (20 tests covering single-level, multi-level, and stratified clustering)
• Enhanced tests to verify functional impact of clustering on standard errors (tests now compare results with/without clustering to ensure standard errors differ appropriately while point estimates remain identical)
• Updated documentation with clustering examples
• Added mathematical description of cluster-robust standard errors to methodology vignette in dedicated subfile (vignettes/articles/_cluster-robust-se.qmd) included after "Finding the MLE numerically" section
  • All mathematical symbols properly defined ($V_{\text{robust}}$, $H$, $B$, $C$, $U_c$, $\nabla_\lambda$)
  • Section header placed in parent file following proper document structure
• Added clustering demonstrations to enteric_fever_example.Rmd vignette using actual SEES Pakistan data with executable examples
  • Cross-references to methodology article for mathematical details
  • Direct comparisons between clustered and non-clustered analyses
  • Side-by-side results showing impact on standard errors while point estimates remain identical
  • Demonstrates single-level clustering with the cluster variable
  • Shows clustering combined with stratified analysis (stratified by catchment with cluster adjustment)
  • Note: Multi-level clustering example was removed as the Pakistan SEES data does not show meaningful differences between single-level and multi-level clustering, making it ineffective for demonstration purposes
• Marked man/ folder as linguist-generated in .gitattributes
• All linting issues resolved
• Version incremented to 1.4.0.9004

Examples

Using with est_seroincidence()

library(dplyr)

# Standard analysis
est1 <- est_seroincidence(
  pop_data = xs_data,
  sr_params = curve,
  noise_params = noise,
  antigen_isos = c("HlyE_IgG", "HlyE_IgA")
)

# With single-level clustered sampling adjustment
est2 <- est_seroincidence(
  pop_data = xs_data,
  sr_params = curve,
  noise_params = noise,
  antigen_isos = c("HlyE_IgG", "HlyE_IgA"),
  cluster_var = "cluster"
)

# With multi-level clustering
est3 <- est_seroincidence(
  pop_data = xs_data,
  sr_params = curve,
  noise_params = noise,
  antigen_isos = c("HlyE_IgG", "HlyE_IgA"),
  cluster_var = c("catchment", "cluster")
)

# Compare results
summary(est1)  # se_type = "standard"
summary(est2)  # se_type = "cluster-robust"
summary(est3)  # se_type = "cluster-robust" (multi-level)

Using with est_seroincidence_by()

# Stratified analysis with clustering
est <- est_seroincidence_by(
  strata = "catchment",
  pop_data = xs_data,
  sr_params = curve,
  noise_params = noise,
  antigen_isos = c("HlyE_IgG", "HlyE_IgA"),
  cluster_var = "cluster"
)

summary(est)
# A tibble: 2 × 14
#   Stratum  catchment     n incidence.rate     SE se_type        
# 1 Stratum… aku          53          0.140 0.0285 cluster-robust
# 2 Stratum… kgh          47          0.200 0.0187 cluster-robust

Point estimates remain unchanged; standard errors appropriately increase to reflect within-cluster correlation (typically 5-15% larger). The se_type column clearly indicates which type of standard error is being used.

Documentation

The methodology vignette now includes a comprehensive section on cluster-robust standard errors (located in vignettes/articles/_cluster-robust-se.qmd and included after the "Finding the MLE numerically" section) explaining:

• Why clustering matters (within-cluster correlation leads to anti-conservative standard errors)
• Mathematical description of the sandwich estimator: $V_{\text{robust}} = H^{-1} B H^{-1}$ with all symbols properly defined
• Implementation examples for single-level and multi-level clustering
• Expected impact on results (point estimates unchanged, standard errors increase 5-15%)

The enteric_fever_example.Rmd vignette now includes an executable section demonstrating:

• How to specify the cluster_var parameter using real SEES data from Pakistan
• Single-level clustering with the cluster variable and comparison with non-clustered analysis
• Combining clustering with stratified analysis (stratified by catchment with clustering adjustment)
• Direct comparisons showing impact on standard errors while point estimates remain identical
• Cross-references to methodology article for mathematical details
• Interpretation of results with actual output
• All examples properly filter both pop_data and noise_params to Pakistan to ensure parameter alignment

Known Limitations

• Parallel processing with clustering in est_seroincidence_by() requires additional work to properly export cluster variables to worker processes. Single-core processing works correctly. A test for this functionality is skipped pending further investigation.

Note on Scope

This PR focuses solely on cluster-robust standard error estimation. Intraclass Correlation Coefficient (ICC) calculation functionality (compute_icc()) that was initially developed has been removed and will be submitted in a separate pull request to maintain a focused scope.

Checklist

• The title of your PR should briefly describe the change.
• Commit/merge messages to be included in NEWS.md should begin with -.
• Code should follow the tidyverse style guide.
• Documentation should use roxygen2, with Markdown syntax.
• Contributions should include unit tests (using testthat).

• Fixes adjusting the standard error for clustered sampling by household/school/etc #471

💡 You can make Copilot smarter by setting up custom instructions, customizing its development environment and configuring Model Context Protocol (MCP) servers. Learn more Copilot coding agent tips in the docs.